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A two-patch prey-predator model with dispersal in predators driven by the 

strength of predation 

Yun Kang^, Sourav Kumar SasmaP, and Komi Messan^ 


Abstract 

Foraging movements of predator play an important role in population dynamics of prey-predator interactions, 
which have been considered as mechanisms that contribute to spatial self-organization of prey and predator. 
In nature, there are many examples of prey-predator interactions where prey is immobile while predator 
disperses between patches non-randomly through different factors such as stimuli following the encounter of 
a prey. In this work, we formulate a Rosenzweig-MacArthur prey-predator two patch model with mobility 
only in predator and the assumption that predators move towards patches with more concentrated prey- 
predator interactions. We provide completed local and global analysis of our model. Our analytical results 
combined with bifurcation diagrams suggest that: (1) dispersal may stabilize or destabilize the coupled 
system; (2) dispersal may generate multiple interior equilibria that lead to rich bistable dynamics or may 
destroy interior equilibria that lead to the extinction of predator in one patch or both patches; (3) Under 
certain conditions, the large dispersal can promote the permanence of the system. In addition, we compare 
the dynamics of our model to the classic two patch model to obtain a better understanding how different 
dispersal strategies may have different impacts on the dynamics and spatial patterns. 

Keywords: Rosenzweig-MacArthur Prey-Predator Model; Self-organization Effects; Dispersal; Persistence; 
Non-random Foraging Movements 


1. Introduction 

Spatial heterogeneity, dispersal patterns, and biotic interactions influence the distribution of species 
within a landscape [4, 46, 59, 67, 68]. Spatial self-organization results from local interactions between organ¬ 
isms and the environment and emerges at patch-scales [60, 73]. For example, limited dispersal ability and 
its related dispersal patterns [56] is considered to be one of the key factors that promotes the development 
of self-organized spatial patterns [I, 41, 69, 70]. 

In nature, especially for ecological communities of insects, dispersal of a predator is usually driven by 
its non-random foraging behavior which can often response to prey-contact stimuli [22], including spatial 
variation in prey density [39] and different type of signals arising directly from prey [74]. For instances, 
bloodsucking insects respond to the carbon dioxide output and the visual signals of a moving animal, which 
in tsetse flies {Glossina spp.) lead to the formation of a “following swarm” associated with herds of grazing 
ungulates [13, 14]. Most mosquitoes were attracted over a larger distance by the odor of the host [16, 17, 18]. 
The wood-wasp, Sirex noctilio, is attracted by the concentration of the scent [18, 52]. Social ants excite 
“pheromone trails” to encourage other individuals to visit the same food source [5]. Plant-feeding insects 
commonly detect food items by gustatory signals [64, 65, 66]. These non-random foraging behaviors driven 
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by prey-mediated patch attractants, prey attractants themselves, and arrestant stimuli following the en¬ 
counter of a prey, can lead to predation rates that are greater in regions where prey are more abundant (i.e., 
density-dependent predation), thus regulate population dynamics of both prey and predator. 

Recent experimental work on population dynamics of immobile Aphids and Coccinellids by Kummel 
et al. [43] show that the foraging movements of predator Coccinellids are combinations of passive diffusion, 
conspecific attraction, and retention on plants with high aphid numbers which is highly dependent on the 
strength of prey-predator interaction. Their study also demonstrates that predation by coccinellids was re¬ 
sponsible for self-organization of aphid colonies. Many ecological systems exhibit similar foraging movements 
of predator. For example, Japanese beetles are attracted to feeding induced plant volatiles and congregate 
where feeding is taking place [51]. Motivated by these field studies, we propose a two-patch prey-predator 
model incorporating foraging movements of predator driven by the strength of prey-predator interaction, to 
explore how this non-random dispersal behavior of predator affect population dynamics of prey and predator. 

Dispersal of predator plays an important role in regulating, stabilizing, or destabilizing population dy¬ 
namics of both prey and predator. There are fair amount literature on mathematical models of prey-predator 
interactions in patchy environments. For example, see work of [46, 24, 3, 23, 25, 26, 63, 72, 19, 20, 61, 37, 38] 
and also see [40] for literature review. Many studies examine how the interactions between patches affect 
the synchronicity of the oscillations in each patch, e.g. see the work of [27, 34], and how interactions may 
stabilize or destabilize the dynamics. For instances, [36, 34] studied a model with two patches, each with 
the well known prey-predator Rosenzweig-McArthur dynamics, linked by density independent dispersal (i.e., 
dispersal is driven by the difference of species’ population densities in two patches). His study showed that 
this type of spatial predator-prey interactions might exhibit self-organization capable of producing stabilizing 
heterogeneities in prey distribution, and spatial populations can be regulated through the interplay of local 
dynamics and migration. 

However, due to the intricacies that arise in density-dependent dispersal models, there are relatively lim¬ 
ited work on models with non-random foraging behavior of predator or non-linear dispersal behavior [39] but 
see the two patch model with predator attraction to prey, e.g. [31]), or predator attraction to conspecific, 
e.g. [57], or only predators migrate who are attracted to regions with concentrated food resources, see the 
work of [21, 7]. Kareiva and Odell [39] proposed a non random foraging PDF model through a mechanistic 
approach to demonstrate that area-restricted search does yield predator aggregation, and explore the the con¬ 
sequences of area-restricted search for predator-prey dynamics. In addition, they provided many supporting 
ecological examples (e.g. Coccinellids, blackbirds, etc.) that abide by their theory. [31] studied a two-patch 
predator-prey Rosenzweig-MacArthur model with nonlinear density-dependent migration in the predator. 
The migration term of the predator is derived by extending the Holling time budget argument to migration. 
Their study showed that the extension of the Holling time budget argument to movement has essential effects 
on the dynamics. By extending the model of [31], Ghosh and Bhattacharyya [15] formulated a similar two 
patch prey-predator model with density-independent migration in prey and density-dependent migration in 
the predator. Their study shows that several foraging parameters such as handling time, dispersal rate can 
have important consequences in stability of prey-predator system. Cressman and Kfivan [9] investigated the 
population-dispersal dynamics for predator-prey interactions in a two patch environment with assumptions 
that both predators and their prey are mobile and their dispersal between patches is directed to the higher 
fitness patch. They proved that such dispersal, irrespectively of its speed, cannot destabilize a locally stable 
predator- prey population equilibrium that corresponds to no movement at all. 

In this paper, we formulates a new version of Rosenzweig-MacArthur two patch prey predator model with 
mobility only in predator. Our model is distinct from others as we assume that the non-random foraging 
movements of predator is driven by the strength of prey-predator interactions, i.e., predators move towards 
patches with more concentrated prey and predator. Our model can apply to many insects systems such as 
Aphids and Coccinellids, Japanese beetles and their host plants, etc. The main focus of our study of such 
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prey-predator interactions in heterogeneous environments is to explore the following ecological questions: 


1. How does our proposed nonlinear density-dependent dispersal of predator stabilize or destabilize the 
system? 

2. How does dispersal of predator affect the extinction and persistence of prey and predator in both 
patches? 

3. How may dispersal promote the coexistence of prey and predator when predator goes extinct in the 
single patch? 

4. What are potential spatial patterns of prey and predator? 

5. How are the effects of our proposed nonlinear density-dependent dispersal of predator on population 
dynamics different from the effects of traditional density-independent dispersal? 

The rest of the paper is organized as follows: In Section (2), we derive our two patch prey-predator model 
and provide its basic dynamics properties. In Section (3), we perform completed local and global dynamics 
of our model, and derive sufficient conditions that lead to the persistence and extinction of predator as well 
as permanence of the model. In Section (4), we perform bifurcation simulations to explore the dynamical 
patterns and compare the dynamics of our model to the traditional model studied by Jansen [35]. In Section 
(5), we conclude our study and discuss the potential future study. The detailed proofs of our analytical 
results are provided in the last section. 

2. Model derivations 

Let Xi, Ui be the population of prey, predator at Patch i, respectively. We consider the following two-patch 
prey-predator interaction model after rescaling (see similar rescaling approaches by [50])(1) 
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where Ki is the relative carrying capacity of prey in the absence of predation; is the relative predation rate 
at Patch 1; di is the relative death rate of predator at Patch i; pi is the relative dispersal rate of predator 
at Patch i; and r is the relative maximum growth rate of prey at Patch 2. All parameters are nonnegative. 
The ecological assumptions of Model (1) can be stated as follows: 

1. In the absence of dispersal, Model (1) is reduced to the following uncoupled Rosenzweig-MacArthur 
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prey-predator single patch models 
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where ri = 1 and r 2 = r and its ecological assumptions [62] can be stated as follows: 

(a) In the absence of predation, population of prey Xi follows the logistic growth model. 

(b) Predator is specialist (i.e., predator yi goes extinct without prey Xi) and the functional response 
between prey and predator follows Hollying Type II functional response. 

2. There is no dispersal in prey species. This assumption fits in many prey-predator (or plant-insects) 
interactions in ecosystems such as Aphid and Ladybugs, Japanese beetles and its feeding plants, etc. 


3. 


The dispersal of predator from Patch i to Patch j is driven by prey-predation interaction strength in 
Patch j termed as attraction strength. In Model (I), we assume that predator in Patch i disperse to 
Patch j is determined by the predation term in Patch j, thus the dispersal term of predator 

from Patch i to Patch j is described by Pi°d^jd^yi which gives the net dispersal of predator at Patch i 
as 


Pi 


1 -I- X, 


Vi 


ajXjyi 
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Vj 


by assuming that the dispersal constant pi of Patch i is the same for predator arriving Patch i from 
other patches as predator leaving Patch i to other patches. This assumption is motivated by the fact 
that dispersal of a predator is usually driven by its non-random foraging behavior which can often 
response to prey-contact stimuli [22] which has been supported in many field studies including the 
recent work by Kummel et al. [43]. 


The state space of Model (1) is K^. Let pi = , vi = ^ 

following theorem regarding the dynamics properties of Model (1): 

Theorem 2.1. Assume all parameters are nonnegative and r,ai, Ki,di,i =1,2 are strictly positive. 
Model (1) is positively invariant and hounded in with lim sup(_,,gQ Xi (t) < Ki for both i = 1,2. In 
addition, it has the following properties: 


1. If there is no dispersal in predator, i.e., Pi = P 2 = 0; then Model (1) is reduced to Model (2) whose 
dynamics can he classified in the following three cases: 


(a) Model (2) always has the extinction equilibrium (0,0) which is a saddle. 

(h) If Pi > Ki or Pi < 0, then the boundary equilibrium {Ki,0) is globally asymptotically stable. 

(c) If ^^ 2 ^ < Pi < Ki, then the boundary equilibrium {Ki, 0) is a saddle while the interior equilibrium 
{pi,Ui) is globally asymptotically stable. 


(d) If 0 < p, < 


Ki-l 


then the boundary equilibrium {Ki,0) is a saddle; the interior equilibrium 


{pi, Vi) is a source, and the system has a unique limit cycle which is globally asymptotically stable. 
In addition, the Hopf bifurcation occurs at pi = . 


2. The sets {(xi, yi, X 2 , j/ 2 ) G : Xi = 0} and {(xi, yi, X 2 , j/ 2 ) G K+ ■ Vi = 0} are invariant for both 

i = 1,2. If Xj = 0, Model (1) is reduced to the single patch model Model (2). If pj = 0, Model (1) is 
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reduced to the following two uncoupled models: 
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where limt_>oo Xj(i) = Kj and the dynamies ofxi,yi is the same as Model (2). 


( 3 ) 


Notes: Theorem 2.1 provides a foundation on our further study of local stability and global dynamics of 
Model (1). In addition, Item 2 of Theorem 2.1 implies that Model (1) has the same the invariant sets Xi = 0 
and Pi = 0 for both i = 1, 2 as the single patch models (2). In addition, the results of the single patch models 
(2) indicate that prey is always persist while predator i is persist ii 0 < pi < Ki hold. 


3 . Mathematical analysis 

Now we start with the boundary equilibria of Model (I). Recall that 

di {Ki - pi){l + pi) r{K 2 - P2){1 + P 2 ) 

ai - di aiKi 02^2 

We define the following notations for all possible boundary equilibria of Model (1): 

^^0000 = (0,0,0,0), E'iciooo = 0,0,0), = {pi, 1 ^ 1 , 0 ,0), £’^^^ 0 ^ 21^2 = 0^ M 2 , ^^ 2 ) 

EK 1 OK 2 O = {Ki,0, K2,0), £’ooif20 = (0,0, £ 2 ,0), = {0,0, p2,i'2), = (mi, M,0) 


The following theorem provides sufficient conditions on the existence and stability of these boundary equi¬ 
libria: 


Theorem 3 . 1 . [Boundary equilibria of Model {l)]Model ( 1 ) always has the following four boundary 
equilibria 

£'0000, £ifi000, £ 00 ^: 20 , EX 1 OK 2 O 

where the first three ones are saddles while EK1OK2O *5 loeally asymptotically stable if pi > Ki and it is a 
saddle if {pi — Ki) {p2 — K2) < 0 or pi < Ki, i = 1 , 2 . Let i,j = 1 , 2 , i j, and 

Ell ~ £/2i!2i 00, £12 “ £/iil2iif20, £21 “ £oO/X2‘^2 0,nd £22 — £^10^12122 ■ 

Then if 0 <pi < Ki, then Model ( 1 ) has additional two boundary equilibria E\i and E’^2 where E[i is always 
a saddle. The boundary equilibrium E[2 is locally asymptotically stable if < pi < Ki and one of the 

following eonditions holds: 


sa: Qj < di,Kj < pj. 

sb: 0 <£, <min{/r„^^}. 


sc ‘ 0 K ■ <r II ■ and o ■ <C ——— ^2 (^3 — ) — 

^ ^ ai-di ^ ^ Mj P] ^ Ui[Kpaj-di)-di]- 

orj• n 11 ■ K _ — and n- ^ _ ^2 1 fl 

sa. V <. pj ^ <. ana pj > , 2 i[di-Kj{aj-di)]- 


And E[2 is a saddle if 0 < pi < ^'2 ^ of the following conditions holds: 
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Kj > max ■ 





0 < ^ < if, 

uc: p,j < Kj < jj-^. 


and pj 


In addition, if 0 < pi < Ki for both i = 1 and i = 2, the boundary equilibria E 12 and E 22 exist but they 
cannot be locally asymptotically stable at the same time while if ri = l,ai = a,di = d,Ki = d for both 
i = 1,2, the boundary equilibria E \2 and E 22 can not be locally asymptotically stable at all if they exist. 


Notes: Theorem 3.1 implies the following points regarding the effects of dispersal in predators: 

1. Dispersal has no effects on the local stability of the boundary equilibrium ExiOK 20 - 

2. Large dispersal of predator in its own patch may have stabilizing effects from the results of Item sd: In 

the absence of dispersal, the dynamics of Patch j is unstable at (if,, 0) since 0 < pj < Kj. However, in 
the presence of dispersal, large values of p, can lead to the local stability of the boundary equilibrium 
E ^2 where i,j = 1, 2 and i ^ j, under conditions of pj < Kj < ■ 

3. Large dispersal of predator in its own patch may have destabilizing effects from the results of Item ub: 

In the absence of dispersal, the dynamics of Patch j is local stable at (if, , 0) since if, < pj. However, 
in the presence of dispersal, large values of p,- can drive the boundary equilibrium Ej 2 being unstable, 
under conditions of 0 < < if, < p, - 

4. Under conditions of pi < Ki, the boundary equilibria £’(’2 and £22 can not be asymptotically stable at 
the same time. 

3.1. Global dynamics 

In this subsection, we focus on the extinction and persistence dynamics of prey and predator of Model 
(I). First we show the following theorem regarding the boundary equilibrium EkiOK^G 

Theorem 3 . 2 . Model (I) has global stability at EK 1 OK 2 O if l^'i > Ki,i = 

Notes: Theorem 3.2 implies that the dispersal of predators does not effect the global stability of the bound¬ 
ary equilibrium EkiOK 20 - 

To proceed the statement and proof of our results on persistence, we provide the definition of persistence 
and permanence as follows: 

Definition 3.1 (Persistence of single species). We say species z is persistent in for Model (1) if 
there exists constants 0 < b < B, such that for any initial condition with z(0) > 0, the following inequality 
holds 

b < liminfz(T) < limsup a;(r) < B. 

r—J-CXD T-j-oo 

where 2 : can be Xi, yi, i = 1,2 for Model (1). 

Definition 3.2 (Permanence of a system). We say Model (1) is permanent in if there exists constants 
0 < b < B, such that for any initial condition taken in with a;i(0)pi(0)a;2(0)y2(0) > 0, the following 
inequality holds 

b < liminf min{a:i(T), j/i(T),a; 2 (T),?/ 2 ('r)} < limsupmax{xi(r), ?/i(t), a; 2 (r), j/ 2 ('r)} < B. 

T—^OC 
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The permanence of Model (1) indicates that all species in the system are persistence. 

Theorem 3.3. [Persistence of prey and predator] Prey Xi,i = 1,2 of Model (1) are always persistent for 
all r > 0. Predator yj is persistent if one of the following inequalities hold 


1. pLj < Kj, pLi > Ki. Or 

2. < Pi < Ki and Kj > max | pj , ^ |. Or 


3. 


Ki-l 


<C 11 ■ K ii ■ <r K ■ <r —^— and o <r — ^ — it — 
^ Pi ^ rvj, pj ana pj ^ „-[di-Kpaj-di)]- 


Notes: Theorem 3.3 indicates that the dispersal of predators does not affect the persistence of preys, while 
small dispersal of predator j, under condition of < pi < Ki, pj < Kj < can keep the persistence 

of predator j. This is consistent with the results of Item uc in Theorem 3.1. 


Theorem 3.4. [Permanence of the two patch dispersal model[ Model (1) is permanent if one of the 
following inequalities hold 


1. < pj < Kj, 0 < <K,<p, and p, > where i = l,j = 2ori = 2,j = 1. 

Or 

2. ^ < h'j < Up Pi P ^' 2 ^ ^ max |/ii, ^ I for both i = l,j = 2 and i = 2, j = 1. Or 


di—Ki (ai—di ) 


3. 


Ki-l 


< Pi < K„ K, > max{^„ < pj < Kj < and Pj < where 


* = l,i = 2 ori = 2,j = 1. 


Notes: According to Theorem 3.3, we can conclude that Model (1) is permanent whenever both predators 
are persistent. Theorem 3.4 provides such sufficient conditions that can guarantee the coexistence of bother 
predator for the two patch model (1), thus provides sufficient conditions of its permanence. Item 1 of this 
theorem implies that if predator j is persistent, and large dispersal of predator i can promote its persistence, 
thus, promote the permanence since in the absence of dispersals in predator, predator i goes extinct due to 
Pi > Ki. This is consistent with the results of Item ub of Theorem 3.1 that large dispersal of predator i can 
have destabilize effects on the boundary equilibria £’^’ 2 . 


3.2. Interior equilibrium and stability 

Let Pi{x) = and qi{x) = ^ have 
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If {x\,y\,X 2 ,y 2 ) is an interior equilibrium of Model (I), then it satisfies the following equations: 


qi{xi) -yi = 0^ qi{xi) = Pi, 

Pi{xi) -di+ pipj {pi{xi) - Pj{xj)) = 0 pi{xi) = 


Piqj{xj)pj(xj)+di 

i+Piqj(xj) 


(4) 


which gives: 


Piqj{xj)pj{xj) + di 


Piqj{xj)pj{xj 


ai(l + Piqjixj)) - {piqj{xj)pj{xj) + df) PiPjixj) [oi - Pj{xj)] Pa^-di 


(5) 
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Since limsup(_j,oc a^i(i) < Ki for both i = 1,2 and yi = qi{xi), therefore, positive solutions of Xi G (0,Ki) 
for (5) determine interior equilibrium of Model (1). By substituting the explicit forms of Pi,qi into (5), we 
obtain the following null dines: 


Xi 


_ a‘2[r‘2plX2(K2-X2)+K2.di] _ 

r2p-LX2{K2ai — K2a2—ax)—r2Pix^{ax—a2)+K2{axr2Pi+aia2—a2di) 




X2 


_ ailriP2Xi{Ki-Xi)+Kid2] _ 

rip2Xi(Kia2 — Kiai—a2)—rip2x‘l(a2—ai)+Ki{a2rip2+aia2—aid2) 


gtjxi) 

gb(xi) 


G(xi) 


( 6 ) 


with ri = l,r 2 = r and the following properties: 

1 . F(0) = = jri --TT = -- T- and FfKz) = 

^ '' fb{0) K2{air2pi-\-aia2-a2di) air’2pi+aia2-a2ai '' fb{f<2) ai-di 

2- ft{x 2 ) = 02 [r 2 PiX 2 {K 2 - X 2 ) + K 2 di] > a 2 K 2 di > 0 for X 2 G [0, K 2 ] and 


fb{x2)\ 


ai—a2—(i 


= o [r 2 Pi{K 2 - X 2 ) + K 2 {a- di)]. 


3. G(0) = ^ = -sw =- ^ and G{Ki) = = M 2 . 

' ' 36 ( 0 ; «i(a2rip2+aia2—ait/2) a2rip2+aia2—01^2 ' ' gb{i<l) 02—02 

4. gt{xi) = Oi [riP 2 a^i {Ki — 2 : 1 ) + Kid 2 ] > aiKidi > 0 for xi G [0, Ki] and 


_ gtjKl) _ d 2 _ 


9bixi)\ 


di —a2 — 


= a [rip 2 {Ki - xi) + Ki{a- ^ 2 )] ■ 


Theorem 3.5. [Interior equilibriumj If pi > Ki for both i = 1,2, then Model (1) has no interior 
equilibrium. Moreover, we have the following two cases: 

1. Assume that Oi > Oj where i = 1, j = 2 or i = 2, j = 1. Define 


Ki (riPj +ai - dj - y^(aj - dj){riPj +0 — dj)^ 


Model (1) has no interior equilibrium if Oi > aj, pi < 


nPj 

AKj aj {ai — aj) {di — ai) 


hold. 


rj(Kjai-Kjaj+ai) 

And it has at least one interior equilibrium (x*, y*, 2 : 2 , j/ 2 ) */ ^^0 following conditions hold for both 
i = 1, j = 2 and i = 2,j = 1 


Oi > max{oj, di, ^ 2 }, Oj > max{di,d 2 }, pj < 


4 :Kiai{aj - ai){dj - Oj) 


F{x 2 ) < Ki, and G{x[) < K 2 


ri {KiOj - KiOi + Qj) 

where sufficient conditions for the inequalities F(xf) < Ki, and G{x'i) < K 2 hold are 

'^{KiUi - Kjdj - dj) 4:Kjaj{Kiai - K^dj - dj) 

~ ^ ^ ajrjK^ + rjKi{Kjaj — KjOi — UiY 




In addition, we have 


^ r- < X* < Ki for both i = 1 ,7 = 2 and i = 2 ,7 = 1. 

-a„-r7.a— 1 J d ^ •d ’ 


2. Assume that ai = 02 = a. Model (1) has no interior equilibrium if di > a + r 2 Pi or d 2 > o + rip 2 
while it has at least one interior equilibrium (a:*, J/*, 3 ^ 212 / 2 ) ^^6 following inequalities hold 


Oi = 02 = o > max{di,d 2 }, F(xf) < Ki, and G(x^) < K 2 





















where sufficient conditions for the inequalities F^x^) < Ki, and G{x1) < K 2 hold are 


Pi < —— 


K^di 

KjVj 


di) 


for both i = l,j = 2 and i = 2, j = 1. In addition, we have ^ ^ < x* < Kj for both i = l,j = 2 

and i = 2, j = 1. 


Notes: Theorem 3.5 provides sufficient conditions on the existence of no interior equilibrium when > Ki 
for either i = 1 or i = 2; and at least one interior equilibrium of Model (1) when < Ki for both i = 1,2. 
The results indicate follows: 

1. It Pi > Ki, then Model (1) has no interior equilibrium if the dispersal of its predator is too small. 

2. It Pi < Ki for both i = 1,2, then large values of the predation rate Oi, aj and small values of dispersal 
of both predators can lead to at least one interior equilibrium. 


The question is how we can solve the explicit form of an interior equilibrium of Model (1). The following 
theorem provides us an example of such interior equilibrium of Model (1). 

Theorem 3.6. [Interior equilibrium and the stability]Suppose that di = d 2 = d. Let 


P-i = 


ai — d’ 


ni = 


{Ki - pi){l + pi) 

aiKi 


,P-2 = 


02 - d 


, 1^2 = 


r{ K2 - P2)it+P2) 
a2K2 


If 0 < Pi < Ki for both i = 1 and i = 2, then E'' = {pi,vi, p 2 ,V 2 ) is an interior equilibrium of Model (1) 
and its stability can be classified in the following cases: 

1. E'^ is locally asymptotically stable if < pi < Ki hold for both i = 1 and i = 2 while it is unstable 
if the following inequality holds 


pi{Kiai — oi — Kid — d) 
aiKi 


rp2{K2a2 - 02 - K2d - d) ^ ^ 
02df2 


2. Assume that [KiOi — oi — Kffi — d)(K 2 a 2 — 02 — K 2 d — d) < 0 and 


pi{Kiai — oi — Kid — d) 
OiKi 


rp 2 {K 2 a 2 - 02 - K 2 d - d) ^ ^ 
02df2 


If KiOi — Oi — Kid — d > 0 (i.e., pi < ^ ) for either i = t or i = 2, then the large values of pi can 

make if* being locally asymptotically stable, i.e.. 


Pi > max ■ 


—Vj — rjpiPj{Kiai — Kid — Ui — d){Kjaj — Kjd — Uj — d) 

[KiKjUjUjdvdai — d)^) ’ 


Kj (L'.iPj-\-l)(aj —d) (Kj^aj^—K.id — a.1 —d) ^ 

Tj fij uiKi{ai — (Kjaj —Kjd — aj —d) 


If, in addition, Oi = 02 = o, Ki = K 2 = K,r = 1, then Pi = P 2 = P = vi = V 2 = v = 
if* = (p, v, p, v) is the only interior equilibrium for Model (1) which has the same local stability as the interior 
equilibrium (p,iy) for the single patch model (2), i.e., if* is locally asymptotically stable if < p < K 
while it is unstable if p < ■ 

Notes: Theorem 3.6 implies Model (1) has an interior equilibrium if* = [pi,vi, p 2 ,i’ 2 ) if di = ^2 = d 
and 0 < Pi < Ki for both i = 1,2. In addition. Theorem 3.6 indicates that dispersal of predators has 
no effects on the local stability if < pi < Ki for both z = 1,2 or one of the single patch models 
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(2) is unstable and + .’~/^ 2 (-R' 2 a 2 ^^g^ K 2 d d) ^ However, large dispersal of predator at 

Patch i can stabilize the interior equilibrium when its single patch model model is unstable at (/i^, Vi) with 

fii{Kiai—ai — Kid—d) . r^2{K2a2 — a2—K2d—d) ^ 

aiKi a2K2 ^ 


4. Effects of dispersal on dynamics 

From mathematical analysis in the previous sections, we can have the following summary regarding the 
effects of dispersal of predators for Model (1): 

1. Large dispersal of predator at Patch i can stabilize or destabilize the boundary equilibrium of x* = 
Ki,y* = 0,x* = iJij,y*j = Vj depending on additional conditions. 

2. Small dispersal of predator at Patch i may preserve its persistence under certain conditions. On the 
other hand, large dispersal of predator at Patch i may promote its persistence when the other predator 
is already persist even if > Ki. 

3. Dispersal has no effects on the persistence of prey and the number of boundary equilibrium. It has 

also no effects on the local stability of the boundary equilibrium s-iid the symmetric interior 

equilibrium {y, v, /i, v) when it exists. 

4. If di > Oi, then small dispersal of predator at Patch i prevents the interior equilibrium while if 
0 < fii < Ki, large predations rates a^, aj and small dispersal of predators at both patches can lead to 
at least one interior equilibrium. 

5. If di = dj and 0 < fii < ^ ^ large dispersal of predator at Patch i can stabilize 

the interior equilibrium {yii,Vi, yj,Vj). 

To continue our study, we will perform bifurcations diagrams and simulations to explore the effects on the 
dynamical patterns and compare dynamics of our model (1) to the classical two patch model (7). 


J^.l. Bifurcation diagrams and simulations 

In this subsection, we perform bifurcation diagrams and simulations to obtain additional insights on the 
effects of dispersal on the dynamics of our proposed two patch model (1). We fix ri = l,r 2 = 1.5,= 
5, K 2 = 3, c?i = 0.2, d 2 = 0.1. Then according to Theorem 2.1, we know that in the absence of dispersal, the 
dynamics of Patch 1 has global stability at (5,0) if 0 < Oi < 0.24; it has global stability at its unique interior 
0.2 (^~ ai°-o.2) (i+ai-0.2) 


equilibrium 


ai-0.2 ’ 


if 2 < 


0.2 

ai — 0.2 


< 5 0.24 < ai < 0.3; and it has a unique limit 


cycle if oi > 0.3; while the dynamics of Patch 2 has global stability at (3, 0) if 0 < 02 < 0.133; it has global 


stability its unique interior equilibrium 


0.1 
02 - 0.1 ’ 




if 1 < 


0.1 

a2—0.1 


< 3 0.133 < 02 < 


0.2 while it has a unique limit cycle if 02 > 0.2. Now we consider the following cases: 


1. Choose oi = 0.25 and 02 = 0.15. In the absence of dispersal, the dynamics at both Patch 1 and 2 
have global stability at its unique interior equilibrium (4,1), (2,1.5), respectively. After turning on 
the dispersal, the coupled two patch model can have one interior equilibrium (see the blue regions in 
Figure 1(a)) which can be locally stable (see the blue dots in Figure 1(b)), or be a saddle (see the green 
dots in Figure 1(b)) where the coupled system has fluctuated dynamics; or it can have two interior 
equilibria (see the red regions in Figure 1(a)) which could be saddles and generate bistability between 
fluctuated interior dynamics and the boundary attractor (4,1,3,0) (see the examples of two interior 
saddles of Figure 4(a)); or it can have three interior equilibria (see the black regions in Figure 1(a)) 
which generate multiple interior attractors (see the examples of two saddles, one sink and two sinks, 
one saddle of Figure 1(b)) or it could have no interior equilibrium (see white and yellow regions of 
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(a) pi V.S. p 2 for the number of interior equilibria (b) pi V.S. predator population yi when p 2 = 

0.025 


Figure 1: One and two bifurcation diagrams of Model (1) where r = 1.5, di = 0.2, d 2 = 0.1, K\ = 5, K 2 = 3, ai = 0.25, and 
a 2 = 0.15. The left figure 1(a) describes how number of interior equilibria changes for different dispersal values pj, i = 1, 2: black 
regions have three interior equilibria; red regions have two interior equilibria; blue regions have unique interior equilibrium; 
yellow regions have no interior equilibrium and predator in Patch 2 dies out; white regions have no interior equilibrium and 
both predator die out. The right figure 1(b) describes the number of interior equilibria and their stability when p 2 = 0.025 and 
pi changes from 0 to 0.5 where y-axis is the population size of predator at Patch 1: Blue represents the sink; green represents 
the saddle; and red represents the source. 


Figure 1(a)). Bifurcation diagrams Figure 1(a), 1(b), and 4(b) suggest that dispersal may destabilize 
system and generate fluctuated dynamics; may generate multiple interior attractors (the case of three 
interior equilibria), thus generate multiple attractors; or even may drive extinction of predator in one 
or both patches (he case of two interior equilibria, no interior equilibrium, respectively). 

2. Choose oi = 0.25 and 02 = 0.25. In the absence of dispersal, the dynamics of Patch 1 has global 
stability at its unique interior equilibrium (4,1) while the dynamics of Patch 2 has a unique stable 
limit cycle around (2,1.5). After turning on the dispersal, the coupled two patch model can have one 
interior equilibrium (see the blue regions in Figure 2(a)) which can be locally stable (see the blue dots 
in Figure 2(b)), or be a saddle (see the green dots in Figure 2(b)), or be a source (see the red dots 
in Figure 2(b)) where the coupled system has fluctuated dynamics for the later two cases; or it can 
have two interior equilibria (see the red regions in Figure 1(a)) which could be two saddles or one 
sink, one saddle and generate bistability between the interior attractor and the boundary attractor 
(see Figure 2(b)); or it can have three interior equilibria (see the black regions in Figure 1(a)) which 
generate multiple interior attractors (see the examples of two sinks and one saddle of Figure 4(b)) or 
it could have no interior equilibrium (see white regions of Figure 2(a)). Bifurcation diagrams of Figure 
2(a), 2(b), and 4(b) suggest that dispersal may stabilize system and generate equilibrium dynamics; 
may generate multiple interior equilibria (the case of three interior equilibria), thus generate multiple 
attractors; or even may drive extinction of predator in one or two both patches (the case of two interior 
equilibria, no interior equilibrium, respectively). 

3. Choose tti = 0.35 and 02 = 0.25. In the absence of dispersal, the dynamics of both Patch 1 and 2 
have a unique stable limit cycle. After turning on the dispersal, the coupled two patch model can 
have one interior equilibrium (see the blue regions in Figure 3(a)) which can be a sink (see the red 
dots in Figure 3(b)) where the coupled system has fluctuated dynamics; or it can have two interior 
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Number of Interior Attractor 



(a) pi V.S. p 2 for the number of interior equilibria (b) p 2 V.S. predator population y\ when pi = 1 


Figure 2: One and two bifurcation diagrams of Model (1) where r = 1.5, d\ = 0.2, d 2 = 0.1, Ki = 5, K 2 = 3, ai = 0.25 and 
a 2 = 0.25. The left figure 2(a) describes how number of interior equilibria changes for different dispersal values pi^i = 1, 2: black 
regions have three interior equilibria; red regions have two interior equilibria; blue regions have unique interior equilibrium; 
yellow regions have no interior equilibrium and predator in Patch 2 dies out; white regions have no interior equilibrium and 
both predator die out. The right figure 2(b) describes the number of interior equilibria and their stability when pi = 1 and p 2 
changes from 0 to 2.5 where y-axis is the population size of predator at Patch 1: Blue represents the sink; green represents the 
saddle; and red represents the source. 


equilibria (see the red regions in Figure 2(a)) which could be two saddles, one sink v.s. one saddle, 
one source v.s. one saddle and generate bistability between the interior attractors and the boundary 
attractor (see Figure 3(b)); it could have no interior equilibrium (see white and yellow regions of Figure 
3(a)). Bifurcation diagrams Figure 2(a)-3(b) suggest that dispersal may generate bistability between 
the interior attractor and the boundary attractor; or even may drive the extinction of predator in one 
or both patches (the case of two interior equilibria, no interior equilibrium, respectively). 

In summary. Figure 1, 2, 3, and 4 suggest that dispersal of predator may stabilize or destabilize interior 
dynamics; it may drive the extinction of predator in one or both patches; and it may generate the following 
patterns of multiple attractors via two or three interior equilibria: 

1. Multiple interior attractors through three interior equilibria: In the presence of dispersal. 

Model (1) can have the following types of interior equilibria and the corresponding dynamics: 

• Two interior sinks and one interior saddle: Depending on the initial conditions with a;i (0)?/i (0)a:2(0)y2(0) > 
0, Model (1) converges to one of two sinks for almost all initial conditions (see examples in Figure 
l(b)-4(b)). 

• One interior sink and two interior saddles: Depending on the initial conditions with a;i(0)?/i(0)a;2(0)j/2(0) > 
0, Model (1) either converges to the sink or has fluctuated dynamics for almost all initial conditions 

(see examples in Figure 1(b), 4(a)). 

We should also expect the case of one sink v.s. one saddle v.s. one source and the case of two source v.s. 
one saddle when the interior sink(s) become unstable and go through Hopf-bifurcation. In addition. 

Model (1) seems to be permanent whenever it processes three interior equilibria. 

2. Boundary attractors and interior attractors through three interior equilibria: 
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Number of Interior Attractor 




(a) pi V.S. p 2 for the number of interior equilibria 


(b) p 2 V.S. predator population y\ when pi = 1 


Figure 3: One and two bifurcation diagrams of Model (1) where r = 1.5, d\ = 0.2, d 2 = 0.1, Ki = 5, K 2 = 3, ai = 0.35 and 
(22 = 0.25. The left figure 3(a) describes how number of interior equilibria changes for different dispersal values pi, i = 1, 2: black 
regions have three interior equilibria; red regions have two interior equilibria; blue regions have unique interior equilibrium; 
yellow regions have no interior equilibrium and predator in Patch 2 dies out; white regions have no interior equilibrium and 
both predator die out. The right figure 3(b) describes the number of interior equilibria and their stability when pi = 1 and p 2 
changes from 0 to 7 where p-axis is the population size of predator at Patch 1: Blue represents the sink; green represents the 
saddle; and red represents the source. 



(a) p 2 V.S. predator population y 2 when pi 
0.5, a2 = 0.15 



P2 


(b) p 2 V.S. predator population y 2 when pi 
0.6, a2 = 0.25 


Figure 4: One dimensional bifurcation diagrams of Model (1) where r = 1.5, di = 0.2, d 2 = 0.1, Ki = 5, K 2 = 3 and ai = 0.25. 
The left figure 4(a) describes describes the number of interior equilibria and their stability when pi = 0.5 and p 2 changes from 
0 to 0.05. The right figure 4(b) describes the number of interior equilibria and their stability when pi = 0.6 and p 2 changes 
from 0 to 1.8. In both figures, blue represents the sink; green represents the saddle; and red represents the source. 
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• one interior sink and one interior saddle: Depending on the initial conditions with xi (0)j/i (0)a:2(0)?/2(0) > 

0, Model (1) converges either to the interior sinks or to the boundary attractors with one predator 
going extinct for almost all initial conditions (see examples in Figure 2(b),3(b), 4(b)). 

• two interior saddles: Depending on the initial conditions with a:i(0)yi(0)x2(0)y2(0) > 0, Model 
(1) converges either to the fluctuated interior attractors or to the boundary attractors with one 
predator going extinct for almost all initial conditions (see examples in Figure 2(b), 3(b), 4). 

• one interior source and one interior saddle: Depending on the initial conditions with xi (0)yi (0)a:2 (0)?/2 (0) > 
0, Model (1) converges either to the fluctuated interior attractors or to the boundary attractors 

with one predator going extinct for almost all initial conditions (see examples in Figure 3(b)). 

Model (1) has bistability between interior attractors and the boundary attractors whenever it processes 
two interior equilibria. This implies that depending on the initial conditions, predator at one patch 
can go extinct when the system has two interior equilibria. 

In general, simulations suggest that Model (1) is permanent when it processes one or three interior equilibria 
while it has bistability between interior attractors and the boundary attractors whenever it processes two 
interior equilibria. 

4-Comparisons to the classic model 

The dispersal of predator in our model is driven by the strength of prey-predator interactions. This is 
different from the classical dispersal model such as Model (7) which has been introduced in [35]: 


dxi 

dt 

dpi 

dt 

dxj 

dt 



l + Xi 


= TjXj 


Ki) l + Xi 
- dipi + pi{pj - Pi) 
L_ 

kJ 1 -k Xj 


dpj 

dt 


1+Xj 


djPj - Pj{Pj - Pi) 


(7) 


where i = l,j = 2 or z = 2,j = 1 with ri = l,r 2 = r. The symmetric case of Model (7) (i.e., = rj, 

Qi = Uj, Ki = Kj, di = dj, and pi = pj) has been discussed and studied by Jansen [35] through simulations 
of different scenarios of local bifurcation analysis. Jansen’s study shows that the classical two-patch model 
(7) has a rich dynamical behavior where spatial predator-prey populations can be regulated through the 
interplay of local dynamics and migration: (i) for very small migration rates the oscillations always synchro¬ 
nize; (ii) For intermediate migration rates the synchronous oscillations are unstable and there are periodic, 
quasi-periodic, and intermittently chaotic attractors with asynchronous dynamics; and (ii) For large predator 
migration rates, attractors in the form of equilibria or limit cycles exist in which one of the patches contains 
no prey. 


Recently, [50] studied Model (7) with both dispersal in prey and predator. Liu provide global stability of 
the interior equilibrium for the symmetric case and performed simulations for the asymmetric cases. Here we 
provide rigorous results on the persistence and permanence conditions that can be used for the comparisons 
to our Model (1) in the following theorem: 

Theorem 4.1. [Summarp of the dpnamics of Model (7)] Define 


Pi = 


di 

[ - di 


Vi = qiifii) = 


ri{K, - i2i)(l + Pi) 
OiKi 


and v) 


Pjdi 
dj + Pj 
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where • Let E\ = = (dij ^^'i, 0, vl) and E 2 = £^ 0 yiX 2 y 2 = (Oj ^1 > d 2 , ^’ 2 ) • d"den we have 

the following summary on the dynamics of Model (7) 

1. Model (7) is positively invariant and bounded in its state space ]R4 

with limsup(_^(^ Xi{t) < Ki for both 

i = 1 and i = 2 . 


2. Boundary equilibria: Model (7) always has the following four boundary equilibria 

Eoooo, EkiOoo, E00K2O, EK1OK2O 

where the first three ones are saddles while EX 1 OK 2 O *5 locally asymptotically stable if 


di + d 2 + Pi + P 2 > 


aiKi 
l + Ki 


02 K 2 

I+K 2 


and 


d aiKx I 

ai - T+i?7 + 


^202 -dl-pi) 

(d2+p2)(l + -ff2) 


> 0 


<;4- 


di — 


o-iKi 


l + Ki 


1 - 


K 2 a >‘2 


{d2+ P2.)(P + K2) 


1 


K2a2 

' d2+P2 

^2 

(1 + K2)\ 


> 0 . 


and it is a saddle if one of the above inequalities does not hold. If 0 < pi < Ki, then the boundary 
equilibrium E^ exists which is locally asymptotically stable if < pt < Ki,rj < Oju'-. 


3. Subsystem i: If xj 

Xi,yi, yj ■' 


0, then Model (7) reduces to the following subsystem (8) with three species 



dyi _ a^Xiyi 

dt 1 + Xi 


Ki) \+Xi 

- d,y^ + p^{yj - yf) 


dVj 

dt 


-djdj - pj{yj - yi) 


( 8 ) 


whose global dynamics can be described as follows: 


3a Prey Xi is persistent for Model (8) with limsup(_^oo Xi(t) < Ki. 

3b Model (8) has global stability at {Ki,0,0) if pi > Ki. 

3c Model (8) has global stability at {pi,Vi,v'^) if < pi < Ki. 

4 . Persistence of prey: Prey Xi persists if Pj < 0, or pj > Kj, or ^ < pj < Kj,ri > OiOj hold. 
Both prey Xi and Xj persist if one of the following three conditions hold 

4(a) Pi > Ki for both i = 1 and i = 2. Or pi < 0 for both i = \ and i = 2. 

4(b) < Pi < Ki,rj > ajVj for both i = l,j = 2 and i = 2,j = 1. 

4(c) Pi > Ki, or Pi < 0, and ^ < pj < Kj,ri > aiHf for either i = l,j = 2 or i = 2, j = 1. 

5. Extinction of prey Prey Xi goes to extinction if ^ < pj < Kj and ^ ^j- dn 

addition, under these conditions. Model (7) has global stability at pj,iyj). 

6. Persistence and extinction of predators: Predator yi and yj have the same persistence and 
extinction conditions. Predators persist if 0 < pi < Ki for both i = \ and i = 2 while both predators go 
extinct if pi > Ki for both i = 1 and i = 2. In addition. Model (7) has global stability at (dCi, 0,7^2, 0) 
for the later case. 
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7. Permanence of Model (7); Model (7) is permanent if 0 < fii < Ki for both i = 1 and i = 2 and 
one of 4(a), 4(b), 4(c) hold. 

8. The symmetric case: Let ri = r 2 = l,ai = 02 = a, bi = b 2 — b, di = d 2 = d and Ki = K 2 = K, 
then yi = fJ .2 = y and vi = V 2 = v. Therefore, we can conclude that Model (8) has global stability at 
{fj., V, y, v) if < y < K. In addition, the loeal stability of {y, v, y, v) for Model (7) is the same as 
the local stability of {y, v) for single patch models when pi = P 2 = 0 0 / Model (7). 

Notes: Theorem 4.1 indicates follows: 

1. li yi > Ki and 0 < yj < Kj, then the large dispersal of predator at Patch i stabilizes EkiQK 20 - 

2. Proper dispersal of predators can drive the extinction of prey in one patch. 

3. Dispersal has no effects on the persistence of predator. This is different from our proposed model (7). 

To see how different types of strategies in dispersal of predators affect population dynamics of prey and 
predator, we start with the comparison of the boundary equilibria of our model (1) and the classic model 
(7). Both Model (1) and (7) always have four boundary equilibria iHooco = (0,0, 0, 0), EkiOoo = {K\,Q, 0,0), 
E 00 K 2 O = (0,0, 7^2,0) and EX 1 OK 2 O = (^ 1 ,0, 7^2,0) among which first three are saddles in both the cases. 
17 0 < yi < Ki for i = 1,2, then Model (1) has another two boundary equilibria E(j^ and E (2 where E(^ is a 
saddle. If 0 < /Xi < Ki, then Model (1) has another boundary equilibrium E(. We summarize and compare 
the dynamics of our model (1) with dispersal in predator driven by the strength of predation and the classical 
model (7) with dispersal in predator driven by the difference of predator densities in Table 1-Table 3. We 
highlight effects of dynamical outcomes due to different dispersal strategies in predators between Model (1) 
and (7) as follows: 

1. The boundary equilibria: EK 1 OK 2 O, -£'i 2 E^. The comparisons listed in Table 1 suggest that 

dispersal of predator has larger effects on the boundary equilibrium of the classic model than ours. 

2. Persistence and extinction of prey. According to the comparison of sufficient conditions leading either 
persistence or extinction of prey in a patch listed in Table 2, we can conclude that the strength of 
dispersal ability of predator has huge impact on the prey for the classical model (7) but not for our 
model (1). 

3. Persistence and extinction of predator. Simulations and the comparison of sufficient conditions leading 
either persistence or extinction of predators in a patch listed in Table 3, suggest that the strength of 
dispersal ability of predator has profound impacts on the persistence of predator for our model (1) 
while it has no effects on the persistence of predator for the classical model (7). 

4. Permanence of a system depends on the persistence of each species involved in the system. Our 
comparisons of sufficient conditions leading to the persistence of prey and predator listed in Table 2-3, 
indicate that dispersal of predator has important impacts in the persistence of predator in our model 
(1) while it has significant effects on the persistence of prey of the classical model (7). We can include 
that (i) the large dispersal of predator in a patch has potential lead to the extinction of prey (the 
classical model (7)) or predator (our model (7)) in that patch, thus destroy the permanence of the 
system; (ii) the small dispersal of predator in Patch i with the large dispersal in Patch j can promote 
the persistence of prey (the classical model (7)) or predator (our model (7)) in Patch i, thus promote 
the permanence of the system. 

5. Interior equilibria: Both our model (1) and the classical model (7) have the maximum number of 
three interior equilibria. However, for the symmetric case, our model (1) can have the unique interior 
equilibrium (see Theorem 3.6) while the classical model can potentially process three interior equilibria 
[35]. 
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Scenarios 

Model (1) whose dispersal is driven by 
the strength of prey-predator interac¬ 
tions 

Classical Model (7) whose dispersal is 
driven by the density of predators 

EK 1 OK 2 O 

LAS and GAS if /i^ > Ki for both 
i = 1,2. Dispersal has no effects on 
its stability. 

GAS if > Ki fo 
LAS if d\ -\- c?2 P 

and l^di j 

Pi Lq 

d 2 +P 2 [ 2 I + K 2 

sal may be able tc 
rium. 

r both 1 = 1,2; While 
+ P2> 

1 02^2 1 

{d,2+P2)ii + K2) \ 

> 0. Large disper- 
stabilize the equilib- 

0 

II 

.og 

LAS if < fii < Ki and one of the 

conditions sa, sb, sc, sd in Theorem 
(3.1) holds. Large dispersal has poten¬ 
tial to either stabilize or stabilize the 
equilibrium. 

Does not exists 

(x. = 0) 

Does not exists 

LAS if < 'fii < Ki and rj < 

ajVj. GAS if < Jii < Ki 

and < i^i- Large disper¬ 

sal of predator in Patch i will either 
destroy or destabilize the equilibrium 
while large dispersal of predator in 
Patch j may stabilize the equilibrium. 


Table 1: The comparison of boundary equilibria between Model (1) and Model (7). LAS refers to the local asymptotical 
stability, and GAS refers to the global stability. 


6. Multiple attractors: Both our model (1) and the classical model (7) have two types bi-stability: (a) 
The boundary attractors where one of prey or predator can not sustain and the interior attractors 
where all four species can co-exist; and (b) Two distinct interior attractors. One big difference we 
observed is that for the symmetric case when each single patch model has global stability at its unique 
interior equilibrium, our model (1) can have only one interior attractor while the classical model can 
potentially have two distinct interior attractors. This is due to the fact that Model (1) has unique 
interior equilibrium while Model (7) can potentially process three interior equilibria as we mentioned 
earlier. 


5. Discussion 

The idea of “metapopulation” originated from Levins [47] where R. Levins used the concept to study 
the dynamics of pests in agricultural field in which insect pests move from site to site through migra¬ 
tions. Since Levin’s work, many mathematical models have been applied to study prey-predator interac¬ 
tions between two or multiples patches that are connected through random dispersion, see examples in 
[36, 35, 6, 42, 34, 45, 46, 58, 8, 25, 2, 54, 28]. The study of these metapopulation models help us get a 
better understanding of the dynamics of species interacting in a heterogeneous environment, and allow us 
to obtain a useful insight of random dispersal effects on the persistence and permanence of these species 
in the ecosystem. Recently, there has been increasing empirical and theoretical work on the non-random 
foraging movements of predators which often responses to prey-contact stimuli such as spatial variation in 
prey density [10, 39], or different type of signals arising directly from prey [74]. See more related examples 
of mathematical models in [48, 44, 11, 7, 12, 21, 9, 43, 15, 31, 27, 55]. Kareiva [40] provided a good review 
on varied mathematical models that deal with dispersal and spatially distributed populations and pointed 
out the needs of including non-random foraging movements in meta-population models. Motivated by this 


17 























Scenarios 

Model (1) whose dispersal is driven by 
the strength of prey-predator interac¬ 
tions 

Classical Model (7) whose dispersal is 
driven by the density of predators 

Persistence of prey 

Always persist, dispersal of predator 
has no effects 

One or both prey persist if conditions 4. 
in Theorem (4.1) holds. Small dispersal 
of predator in Patch i and large disper¬ 
sal of predator in Patch j can help the 
persistence of prey in Patch i. 

Extinction of prey 

Never extinct 

Xi extinct if ^ < 'jlj < Kj and 

^ Large dispersal of 

predator in Patch i can promote the ex¬ 
tinction of prey in Patch i. 


Table 2: The comparison of prey persistence and extinction between Model (1) and Model (7) 


Scenarios 

Model (1) whose dispersal is driven by 
the strength of prey-predator interac¬ 
tions 

Classical Model (7) whose dispersal is 
driven by the density of predators 

Persistence of predator 

Predator at Patch j is persistent if Con¬ 
ditions in Theorem (3.3) holds. Small 
dispersal of predator in Patch j can 
help the persistence of predator in that 
patch. Dispersal is able to promote the 
persistence of predator when predator 
goes extinct in the single patch model. 

Predators in both patches have the 
same persistence conditions. They per¬ 
sist a 0 < Hi < Ki for i = 1,2. Dis¬ 
persal seems to have no effects in the 
persistence of predator. 

Extinction of predator 

Simulations suggestions (see the yellow 
regions of Figure 1(a) and Figure 3(a)) 
that the large dispersal of predator in 
Patch i may lead to the its own extinc¬ 
tion. 

Predators in both patches have the 
same extinction conditions. They go 
extinct if Hi > Ki or Hi < 0 for f = 1, 2. 


Table 3: The comparison of predator persistence and extinction between Model (1) and Model (7). 


and the recent experimental work of immobile Aphids and Coccinellids by [43], we formulate a two patch 
prey-predator model (1) with the following assumptions: (a) In the absence of dispersal the model reduced 
to the two uncoupled Rosenzweig-MacArthur prey-predator single patch models (2); (b) Prey is immobile; 
and (c) Predator foraging movements are driven by the strength of prey-predator interaction. We provide 
basic dynamical properties such positivity and boundedness of our model in Theorem 2.1. 

Based on our analytic results and bifurcation diagrams, we list our main findings regarding the following 
questions stated in the introduction how our proposed nonlinear density-dependent dispersal of predator 
stabilizes or destabilizes the system; how it affects the extinction and persistence of prey and predator in 
both patches; how it may promote the coexistence ; and how it can generate spatial population patterns of 
prey and predator: 

1. Theorem (3.1) provides us the existence and local stability features of the eight boundary equilibria 
of our model (1). This result indicates that large dispersal of predator in its own patch may have 
both stabilizing and destabilizing effects on the boundary equilibrium depending on certain conditions. 
Theorem (3.2) gives sufficient conditions on the extinction of predator in both patches, which suggest 
that predator can not survive in the coupled system if predator is not able to survive at its single 
patch. In this case, dispersal of predator has no effect on promoting the persistence of predator but 
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dispersal may drive predator extinct even if predator is able to persist at the single patch state (see 
white regions of Figure 1(a), 2(a), and 3(a)). 

2. Theorem (3.3) provides sufficient conditions of the persistence of prey and predator while Theorem 
(3.4) provides sufficient conditions of the permanence of our two patch model. These results imply that 
under certain conditions, large dispersal of predator can promote its persistence, thus, promote the 
permanence of the coupled system while predator in that patch goes extinct in the absence of dispersal. 
Our numerical studies also suggests that large dispersal can also drive the extinction of predators in 
both patches (see white regions of Figure 1(a), 2(a), and 3(a)). 

3. Theorem (3.5) and Theorem (3.6) provide sufficient conditions on the existence and the local stability 
of the interior equilibria under certain conditions. Our analytic study shows that large dispersal of 
predator may be able to stabilize the interior equilibrium when one of the single patch has global 
stable interior equilibrium while the other one has limit cycle dynamics. At the mean time, our 
bifurcation diagrams (see Figure 1(b), 2(b), and 2(b)) suggest that the stabilizing or destabilizing 
effects of predator’s dispersal are not definite, i.e., dispersal can either stabilize or destabilize the 
system depending on other life history parameters. Moreover, our simulations also suggest that the 
dispersal of predator can either generate multiple interior equilibria or destroy the interior equilibrium 
which leads to the extinction of predator in one patch or predators in both patches. 

Comparisons to the classic model (7): We provide detailed comparison between the dynamics of our 
model (1) to the classic model (7). These comparisons suggest that the mode of forging movement of preda¬ 
tor has profound impacts on the dynamics of the coupled two patch model. Here we highlight two significant 
differences: (1) the strength of dispersal ability of predator has profound impacts on the persistence of 
predator for our model (1) while it has no effects on the persistence of predator for the classical model (7). 
However, the dispersal of predator has huge impacts on the persistence of prey for the classical model (7) 
while it has little or no effects on the persistence of prey for our model (1). And (2) for the symmetric case, 
our model (1) has a unique interior equilibrium while the classical model (7) can have up to three interior 
equilibria thus it is able to generate different spatial patterns. 


Future work: Our study combined with the literature study on the classical model (7) by [50, 36, 35, 25, 46], 
provide us a better understanding on how different dispersal behavior of predator could have different effects 
on the dynamical outcomes and spatial pattens. In nature, predator may have different foraging behavior as 
pointed out by Kummel et al. [43] that the foraging movements of predator Coccinellids are combinations 
of passive diffusion, conspecific attraction, and retention on plants with high aphid numbers. It will be 
interesting to the extended version of Model (1) and Model (7) by incorporating two different modes of 
foraging behavior. One potential example is showed as follows: 


dxi 

dt 

dyi 

dt 

dX2 

dt 

dy2 

dt 


Xi 



aixiyi 

l+Xi 


aiXiyi 

1 -I- Xi 


diyi -l-pis 


aiXiyi 

1 -l- 2:1 


2/2 


a2X2y2 

1+ X2 



rx2 



a2X2V2 
1 + X2 


a2X2y2 

1+ X2 


^22/2 + /O 2 S 


a2X2y2 

I+X 2 


2/1 


aiXiyi 

l+Xi 



Pl(l-s) ( 1 / 2 - 2 / 1 ) 


P 2 (l-S) (2/1 -2/2) 


(9) 


where s is a real number in [0,1] indicating the portion of predator using the dispersal strategy driven by 
the strength of the predation (our model (1)) and 1 — s indicates the portion of predator using the dispersal 
strategy driven by the density difference of predator in two patches (the classical model (7)). It will be 
even more interesting to develop a two patch model with adaptive dispersal strategies by letting s change 
over time and depend on the fitness of predator. These are ongoing research projects by the authors. We 
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would also like to point out that Bruder et al. (Bruder, A., Thompson, H., Brown, D. & Kummel, M: 
Pattern formation in a two-patch predator prey-model with diffusion and attraction to predation, working 
on progressing) is working on Model (9) focusing on the spatial patterns generated by these two strategies 
of dispersal in predator. 


6. Proofs 


Proof of Theorem 2.1 


Proof Notice that both = 0 and = 0 for j = 1,2, thus according to Theorem A.4 (p.423) 

in [71], we can conclude that the model (1) is positive invariant in Now we can go ahead to show 
the boundedness of the system. First, we have the following inequalities due to the property of positive 
invariance: 


dxi ( Xi 


a^Xiyi ^ 

—— < nx^ 1 - ^ 

IPXi \ Kt 


which implies that 


limsupa;i(t) < Ki. 

t—¥oo 


Define V = P 2 {xi + yi) + pi{x 2 + 2 / 2 ), then we have 


dV 

dt 


_ „ d(a;i-|-yi) , ^ d{x2+y2) 

- P‘^ - di -- di - 

= p2a:i (l - j -h pirx2 f 1 - - P2diyi - pid2y2 

= P2X1 (^1 - j -H Pirx2 (^1 - j -H P2diXi + pid2X2 - P2di{xi + yi) - Pid2{x2 + 2/2) 

< M - d[p2{,xi P yi) P pi{x2 + y2)] = M - dV 


where d = min{(ii,c? 2 } and 

M= max (p 2 a;i ( 1 -+ di I I + niax (pia :2 ( 1 - ^ ^2 | )• 

o<xi<ifi L \ Ki ) ) q<x2<K2 L \ K 2 J J 

Therefore, we have 


lim sup V (t) = lim sup p 2 (xi (t) + 2/1 (t) ) -f pi (x 2 (t) -f 2/2 (f)) < -y 

t—¥oo t—¥oo 

which implies that Model (2) is bounded in 

If there is no dispersal in predator, i.e., pi = 0,i = 1,2, we can easily check that Model (1) is reduced 
to the two uncoupled Rosenzweig-MacArthur prey-predator single patch models (2) with ri = 1 and r 2 = r. 
The global dynamics of the single patch model (2) can be summarized from the work of [49, 50, 29, 30]. 
Thus, we omit the detailed proof here. 

Recall that both = 0 and ^|y._Q = 0 for z = 1,2, therefore, the sets {(xi, 2/1,3:2,2/2) £ K+ : 

Xi = 0} and {(xi, 2 / 1 , X2,2/2) £ IR+ ■ yi = 0} are invariant. This indicates that if Xj{ 0 ) = 0, then Xj{t) = 0 for 
all t > 0. Therefore, the population of yj converges to 0 since 

^ = -djVj - ^ 0 ^ 

Applying the results in [53], we can conclude that Model (1) is reduced to the single patch model (2) when 
Xj = 0. 

In the case that yj = 0, Model (1) is reduced to Model (3) by replacing yj = 0 in Model (1). 
Summarizing the discussions above, we can conclude that the statement of Theorem 2.1 holds. □ 
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Proof of Theorem 3.1 

Proof According Theorem 2.1, sufficient condition for the single patch model (2) having the unique inte¬ 
rior equilibrium {fjLi,Vi),i = 1,2 is /i^ < Ki. Therefore, sufficient condition for Model (1) having boundary 
equilibria is Mi < ^i- Similarly, sufficient condition for Model (1) having boundary 

equilibria £' 00 ^ 21^2 and Eki 0 p. 2 U 2 is M 2 < £ 2 - 

The local stability of an equilibrium (a;^, y*, can be determined by the eigenvalues A^, * = 1, 2, 3, 4 

of the Jacobian matrix J(xi ,yi ,x* ,y*) (10) of Model (1) evaluated at the equilibrium. 






2x1 \ 

ai!/*(l + Pl!/ 2 ) 

(1 + ^ 1 )^ 

0 


_P2fMJCJ^ 

(1 + ^ 0 ^ 


1+x* 


aia;*(l + piBj) 
l + iJ 


pia2X2y2 

1+x'i 


-di 


P 22/2 


( <^2^2 aix^ \ 

I + X 2 l+a^i J 




pio-2yiy2 

(l+“2)^ 

2x. 


<^2^2 


“2!/2(l+P2ai) 

( 1 +^ 2 )^ 




0.2X<2 

1+x* 


021:2 (1 + P2B1) 
1 + 2 ;; 


P20ia:i!/i i 

1 + 2:1 “2 


After substituting the boundary equilibria AqooOj AiCiOOOi £ 00 X 20 ) £aiipi 00 and £ 00 / 221^2 into the Jacobian 
Matrix (10), we can conclude that these equilibria are saddles since they have both positive and negative 
eigenvalues. 


The eigenvalues of (10) evaluated at £xi 0 X 20 are as follows: 


Ai 


-1(<0), A 2 


QiATi 
1-f £1 


dl < 0 Ml ^ -£l) 


-r (< 0 ), A 4 


02 £2 


d 2 <0^ H 2 > £ 2 - 


Therefore, £xi 0 X 20 is locally asymptotically stable if Mi > Ei,i = 1,2 while it is a saddle if either 
(Mi — £ 1 ) (m 2 — £ 2 ) < 0 or < £i, f = 1 , 2 holds. 


Now we focus on the local stability of £;xiPiX 20 and £xi 0 ai 2 P 2 when they exist. After substituting the 
boundary equilibrium £^ 1121 X 20 to (10), we can obtain the eigenvalues of the Jacobian matrix evaluated at 
this boundary equilibrium as follows: 


and 


Ai = -r, A2 = 


K2{a2-d2)-d2 [£ 2(02 - di) - di] 


1 + £2 


M2 


(1 + £ 2 ) 


A3 J- A4 


£i(ai - di) - (oi + di) 
ai£i(ai - di) 


A3A4 


di(£i(ai — di) — di) 

Oi£i 


Notice that the eigenvalues of A3 and A4 being negative is equivalent to the case that the unique interior 
equilibrium (^ 1 , vi) being locally asymptotically stable for the single patch model (2) when i = 1. Thus, we 
can conclude that < mi < £1 are sufficient conditions for A 3 and A 4 being negative. Now we explore 

sufficient conditions for A 2 being negative. First, we have mi < £1 due to the existence of Epi^^p^K 20 - We 
have the following three cases: 


1. If M 2 > £2 ^ £ 2(02 — d 2 ) — d 2 < 0, then the first term of A 2 is negative. This also implies that Model 
(1) has no boundary equilibria of £ 00 / 22^2 and £xi 0 / 22^2 • Since mi < £i) therefore, we have A 2 < 0 for 
all p 2 > 0 if 02 < di or £2 < since 


£ 2(02 — di) — di < 0 either 02 < di or £2 < -—. 

02 - di 


( 10 ) 
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Therefore, we can conclude that A 2 is negative if either 02 < di, K 2 < ^2 or K 2 < min |/i 2 , | ■ 

Assume that ^ 2(02 — di) — di > 0 (i.e., K 2 > > 0), then A 2 < 0 if p 2 is small enough, i.e., 

satisfies the condition of p 2 < case, we can conclude that A 2 is negative if 

d2 - K2{a2 - d2) 


0 < 


di 

02 — di 


< K 2 < P 2 and p 2 < 


vi [K2{a2 - di) - di] ■ 


2. If p 2 < K 2 K 2 {a 2 — d 2 ) — d 2 > 0, then the first term of A 2 is positive. This also implies that Model 
(1) has two boundary equilibria of £ 1 oom 21^2 and £'ifi 0 / 22 i^ 2 - this case, sufficient conditions for A 2 
being negative are K 2 < and p 2 large enough. More specifically, p 2 has to satisfy the following 

inequality: 

1 ^ 2(02 — d2) — d2 
^ 1^1 [di — K2{a2 — di)] 

Therefore, A 2 is negative if p 2 < K 2 < and p 2 > ■ 

Summarizing the discussions above, we can conclude that the boundary equilibrium is locally 

asymptotically stable if < Pi < Ki and one of the following conditions holds: 

1. 02 < di,K2 < P2- 

2. 1^2 <min{/r2,^^}. 

3. 0 < ^ <K 2 <P 2 and p 2 < ■ 

4. P 2 <K 2 <^^ and p 2 > ■ 

And is a saddle if pi < following conditions holds: 


1. 1^2 >max{Ai2,y^}. 

2- 0 < ^ < < M 2 and p 2 > ■ 


3. P 2 <K 2 < ^ and p 2 < ■ 

Similarly, we can obtain sufficient conditions for the local stability of the boundary equilibrium ExiO^i^v^ as 
the statement. 

If Pi < Ki, then Model (1) has the boundary equilibria and EKi 0 fj, 2 i ^2 according to Theorem 3.1 

and the discussions above. If both and EKiOfj .2122 are locally stable, then the following inequalities 

are satisfied: 

^2 di *^2 j j 

Ml < ^ -p < -p ^ di < d2 


and 


M 2 < K 2 < 


oi — d2 oi — di ^ oi — d2 
di d2 ^ di 


d2 < di 


02 — di 02 — d2 02 — di 
which are contradiction. Therefore, E^-^^^^k 2 Q and EKi 0 fj. 2 V 2 can not be local stable at the same time. 

Now if Vi = 1, Oi = o, di = d, Ki = d for both i = 1,2, then E^ 2 ^ * = Ij 2 exist if Ka — Kd — d > 0. This 
implies that one of the eigenvalues of the Jacobian matrix of Model (1) evaluated at E ^2 is positive, i.e.. 


K{a — d)^Ka — Kd — d) + pj{Ka — Kd — d)^ 

A:(A: + l)(o-d)2 


> 0 


which indicates that E ’^2 can not be stable for both i = 1 and i = 2. 


□ 
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Proof of Theorem 3.2 

Proof Let Pi{x) = and qi{x) = ^i(Ki-oc)(i+x) ^ 


r.i. I 1 - I - 


aiX^y^ 


QjjXi 


(l+Xi) l + Xi 


ri{Ki - a;i)(l + Xi) 


= Pi{xi) [qi{xi) - yi\. 


We construct the following Lyapunov functions 


nx 

Viixi.yi) =P2 

JKi 


and 


V 2 {x 2 ,y 2 ) = Pi 

Jh 


PliO -PijKi) 

PiiO 

P2iO-P2{K2) 


df + P2yi 


df + Piy 2 


IK 2 P2(C) 

Now taking derivatives of the functions (11) and (12) with respect to time t, we get 


( 11 ) 

( 12 ) 




Pi{xi) 

P2 \piixi) -pi{Ki)] [qi{xi) - yi] + P2yi [pi{xi) - di] + PlP2yiy2 [pi{xi) - P2{X2)] 
P2 \piixi) -pi{Ki)] qi{xi) + P2yi [pi{Ki) - di] + PlP2yiy2 [Plixi) -P 2 ix 2 )] 

(13) 


and 


d T/ fj.\ _ „ P2{x2)-P2iK2) dX2 i „ dy2 

jjV2[X2(t),y2[t)) - Pi -- dT+P^-dT 

= Pi [P2 {x 2) -P2{K2)] [q2{x2) - 1 / 2 ] + Piy2 [p2{x2) - ^ 2 ] + PiP22/i2/2 [^2( 2 ^ 2 ) “ Pi(a;i)] • 
= Pi [P 2 (a; 2 ) - P 2 {K 2 )\q 2 {x 2 ) + Piy 2 [^ 2 ( 1 ^ 2 ) - ^ 2 ] + PiP 2 yiy 2 \P 2 {x 2 ) -pi(a;i)] 

(14) 

Let V = Vi + ¥ 2 - Now adding (13) and (14), we get 

= ^Vi{xi{f),yi{t)) + ^V 2 {x 2 it),y 2 {t)) = p 2 [piixi) - pi{Ki)] qi{xi) + p 2 yi [piiKi) - di] 

+ Pi [P2ix2) - P2{K2)] q2{x2) + Pl?/2 [P2(1^2) “ ^ 2 ] ■ 


Since pi > Ki a ‘^d > Ki i^ — Pii^i) < di, thus, we have Pi{Ki) — < 0. Notice that Pi{xi) is an 

increasing function in Xi, therefore, \pi{xi) — Pi^Ki)] is positive for Xi > Ki and it is negative for Xi < Ki. 
At the mean time, we have qi{xi) is positive for Xi < Ki and it is negative for Xi > Ki. This implies that 
[Pi{xi) — pi{Ki)\ qi{xi) < 0 for all Xi > 0. Therefore, we have < 0 in This implies that both Vi and 
V 2 are Lyapunov functions, and the boundary equilibrium EK 1 OK 2 O = (ATi, 0, ^ 2 ,0) is globally stable when 
Pi > Ki according to Theorem 3.2 in [30]. □ 


Proof of Theorem 3.3 

Proof. According to Theorem 2.1, we know that Model (1) is attracted to a compact set C in K^. Moreover, 
if yj — 0, Model (1) is reduced to the two uncoupled models (3) while if Xj = 0 it is reduced to a single 
patch model (2). 

First we focus on the persistence conditions for prey xi. Model (1) is reduced to a single patch model 
(2) when a:i(0) = 0, i.e., we have xi = yi = 0. Notice that 

xi \ _ aiyi 
kJ 1+xi 

According to Theorem 2.5 of [32], we can conclude that prey Xi is persistent. Similarly, we can show that 
prey xi is persistent for all r > 0. 


= 1 > 0 . 




dxi 


xidt 


xl—0,yi—0 
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Since both Xi and X 2 are persistent, then we can conclude that Model (1) is attracted to a subcompact 
set Cs of C that excludes i^oooo, E 00 K 2 O and £ 00 X 20 - Therefore, we can restrict the dynamics of Model (1) 
on the compact set Cs- Now we focus on the persistence conditions for predator yi. According to Theorem 
2.1, if 2/1 = 0, Model (1) is reduced to the two uncoupled models (3). In this case, according to both 
Theorem 2.1 and 3.1, the omega limit sets of (1) on the compact set Cs are Ek^oog, , EK 1 OK 2 O, EKi 0 ti 2 i ^2 if 
< /i 2 < K 2 while they are EK 1 OK 2 O if M 2 > K 2 - Now we consider the following two cases: 

1. If 2^2 > K 2 , according to Theorem 2.5 of [32], we can conclude that predator yi is persistent if all of 
the following equations are strictly positive: 


dyi 

yidt 


EkiOOO 


dyi 

Vldt 




aiXi 

l+tci 

aiXi 

1+Xi 


^ \ ^ (aixiy2 

a-i +Pi 

^ \ ^ (0'i^iy2 

di +Pi 


CL2X2y2 \ 


o,K, 

1 + 2:2 ) 


P 1 + ^1 

^K-iQQO 

CL2X2y2 A 


a\K\ 

1 + 2:2 ) 


rp l-\-Ki 

J^KiOK20 


Since — di > 0 pi < Ki, therefore, we can conclude that predator yi is persistent if yii < Ki 

and y ,2 > K 2 - 

2. If < p 2 < K 2 , according to Theorem 2.5 in [32] and discussions above, we can conclude that 

predator yi is persistent if ni < Ki and the following equation is strictly positive: 


dyi 

yidt 




aixi 
1 + Xi 

1+Ki 


— di -\- Pi 


( aixiy2 
\ 1+xi 


-di + P1V2 (ypfy - 


a2X2V2 \ 

1 + 2:2 J 


+'K"iO/^2^2 


0-2y-2 1 _ a-iXi 

1 + /22 ) 1 + Ki 


_ Old^l 

~ 1+Ki 


di + Pi 


( aiKiV2 

V 1+^fl 


di + P1V2 -d-^ > 0 


02112 02 \ 
1 + 1^2 J 


According to the proof of Theorem 3.1, we can see that sufficient condition that 


dyi 

Vldt 


> 0 


E 


K’i0m2^'2 


holds is the same as sufficient condition for the boundary equilibrium EKi 0 y. 2 V 2 being unstable when 
Pi < Ki. Therefore, we can conclude that predator yi is persistent if one of the following inequalities 
hold 


(a) pj < Kj,p^ > Ki. Or 

(b) < Pi < Ki and Kj > max ^pj, Or 

^ < Pi < Ki, Pj < Kj < yrCj. and pj < , 2 i[di-Kj{ai-di )]' 

where i = l,j = 2orz = 2,j = l.one of the following inequalities hold. 


Based on the discussion above, we can conclude that the statement of Theorem 3.3 holds. 

□ 


Proof of Theorem 3.4 

Proof. If < Kj,pi > Ki, then according to Theorem 3.3, we can conclude that prey Xi for both 

i = 1,2 and predator yj is persistent. This implies that Model (1) is permanent if predator yi is persistent. 
Since ExfE < < Kj, then Theorem 3.1 indicates that the omega limit set of Model (1) when = 0 is 

Efi^viK 20 when i = 2, j = 1 while its omega limit set is EKi 0 fi 2 V 2 when i = 2,j = 1. Now let i = 1, j = 2, then 
according to Theorem 2.5 of [32], we can conclude that predator yi is persistent if the following equation is 
strictly positive: 


dyi 

Vldt 


diiKiOiJ.21'2 


f^-dl+PlE2 


a2X2y2 \ 
1 + 2:2 / 

Q 2/22 ) 

1+/22 J 


diiKi0ii2O2 

oiEi _ 

1+ifl 


_ oiEi 

~ 1 + Ki 


di + Pi 


( OlKi V2 

V 1+ifi 


di + P1V2 - ^ 2 ) > 0 


0211202 A 
1 + 02 ) 
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Since fj,i > Ki i+^\ — di <0, therefore, sufficient conditions for 


yidt 


EKiOfJ.21'2 


> 0 is 


aii^i J n ^ ^2 , di-Ki{ai-di) 

— 02 > 0 Ai > - — and pi > 


l + Ki 


Oi — d2 


V2 [Ki{ai - d2) - d2]' 


Similarly, we can show that predator 7/2 is persistent when i = 2,j = 1. Therefore, Model (1) is permanent 
if the following inequalities hold for either i = 2, j = 1 or i = 1, j = 2, 

^7 ~ 1 ^ c?7 , di — KAtti — di) 

-A -— < pj < > Ki,0 < -—and p* > 


Oi — dj 


Vj [Ki{ai - dj) - dj] ■ 


According to Theorem 3.3, we can conclude that prey Xi for both i = 1,2 and predator yi is persistent 
if the following inequalities hold 


K,-l 


^ fdi ^ , Pj 


- 1 .7 

—- and Ki > max 

o J 


o, - d, }• 


dj di 


Therefore, Model (1) is permanent if the above inequalities hold for both i = 1, j = 2 and i = 2,j = 1. On 
the other hand, predator yj is persistent if the following inequalities hold 


K,-l Kj - 1 

< Pi < Ki, 


di j 

< pj < Kj < -— and pj < 

d d (X j 


~ dj) — dj 

Vi [di - Kjidj - di)]' 


Therefore, both predator yi and yj are persistent if the following inequalities hold for either z = 1, j = 2 or 
i = 2,j = 1, 


K,,-l 


< Pi < Ki, Kt > max (/Zj, — 

t di — di J 


1 ^ ^ rv ^ _J „ ^ dAj{dj dj) dj 

<Pj<Kj< and pj < 


dj — di 


Vi [di - Kj{dj - di)]' 


Based on the discussion above, we can conclude that the statement of Theorem 3.4 holds. 


□ 


Proof of Theorem 3.5 

Proof. If Pi > Ki for both z = 1, 2, then Model (1) has global stability at (ATi, 0, ^^ 2 ,0) according to Theorem 
3.2. This implies that Model (1) has no interior equilibrium when pi > Ki for both z = 1, 2. 

The interior equilibrium (xj, j/J, z/2) is determined by the positive intersections of the nullclines (6) 
a;i = F{x 2 ) = and X 2 = G{xi) = where 

ftix2) = d2 [r2PiX2 {K2 - X2) + K2di], fb{x2) = r2PiX2 {K2ai - K2d2 - di)-r2Pixl{di-d2)+K2{dir2PiAdid2-a2di) 


and 


gt{xi) = Oi [rip 2 a;i (A'l - xi) + * 11 ^ 2 ], 9b{xi) = riP 2 Xi {Kid 2 - Kjdi - d 2 )-riP 2 xl{d 2 -di)+Ki{d 2 riP 2 +did 2 -aid 2 ). 
Notice that the nullclines Xi = F{x 2 ) = and X 2 = G{xi) = has the following properties: 


1. F(0) = ^ ^ ^^ and F(K 2 ) = 44^ = = pi- 

i > Jb(0) K2(air2pi+aia2—a2di) air2Pi+aia2 —0201 v Jb{J<2) ai—di 

2. ft{x2) = 02 [r 2 PiX 2 {K2 - X2) + A: 2 di] > d2K2di > 0 for X2 G [0, K2] and 


di 


fb{x 2 )\ 


di —< 22 — O' 


3- G(0) = = 


aiKid2 


pb(0) Ki{a2rip2+aia2—a-id2) a2rip2-\-aia2—aid2 


= d [r 2 Pi{K 2 - X2) + K 2 {d - di)]. 

- “d G(A',) = Sis = = ». 
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4. gt{xi) = ai [rip 2 Xi {Ki — xi) + Kid 2 \ > aiKidi > 0 for xi G [0, Ki] and 

5b(2^l)|ai=a2=a = “ [’■lP2(^l “ a^l) + Ki{a - ^2)] . 


According to Theorem 2.1, we know that population of prey Xi for i = 1,2 has the following properties: 

limsupa;i(t) < Ki. 

t—¥00 

Thus, we can restrict the function F{x 2 ) on the domain of [0, K 2 ] and G{xi) on the domain of [0, Ki], Since 
ft{x2) > 02^12^1 > 0 for X2 G [0,if2] and gt{xi) > aiKidi > 0 for Xi G [0, ATi], thus, Model (1) has no 
interior equilibrium if fb{x 2 ) < 0 for X 2 G [0,^^ 2 ] or gb(xi) < 0 for xi G [0, ATi]. 


Now we assume that oi > 02 , then we have fb{x 2 ) < 0 for X 2 G [0, K 2 ] if 
rlpi {K2ai — K2a2 — ai)^ < AK2r2pi{a2 — ai)(air2pi + 0102 — a2di) pi < 
while fb{x 2 ) > 0 for X 2 G [ 0 ,^ 2 ] if ai > since 


41^202(01 — a 2 ){di — ai) 
r 2 {K 2 ai — K 2 a 2 + oi)^ 


/&(0) = K 2 {air 2 Pi + 02(01 - di)) > 0 and fb{K 2 ) = 02 ^: 2(01 - di) > 0. 


And gb{xi) > 0 for Xi G [0, ATi] if 


riP 2 {Kia 2 - Aioi - 02)^ < 4 Airip 2 (ai - a 2 )(a 2 rip 2 + 0102 - 01^2) P 2 < 


4Aiai(ai — 02) (02 — (^2) 
n {Kia 2 - Kiai + a 2 )^ 


Similar cases can be made for oi <02, therefore we can conclude that Model ( 1 ) has no interior equilib¬ 
rium if either 

4^202(01 - 02)(dl - Oi) 

Ol > 02, Pi < ---^ 

r2 {K 2 ai — K 2 a 2 + oi) 

or 

4AriOi(o2 — oi)(d2 — 02) 

O2 > Ol, P2 < ---^ 

ri {Kia2 - ATiOi -I- 02) 

hold. 

Now we focus on sufficient conditions lead to both fb{x2) > 0 for X2 G [0,^12] and gb{xi) > 0 for 

xi G [ 0 , ATi]. Assume that oi > 02 and Oi > di for both z = 1 , 2 . Notice that /b( 0 ) = K2{air2Pi -I- 02(01 — 

di)) > fb{K2) = 02^12(01 — di) > 0 with 

fb{x2) = r2PiX2 {K2ai - ^202 - oi) - r2Pi0:2(01 - 02) -I- K2{air2Pi + 0102 - 02di). 

Therefore, we have fb{x2) > fb{K2) for all 0:2 G [ 0 , ^12]- Since gb{xi) is a degree 2 polynomial with the 
positive coefficient in the degree 2 and 5&(0) = Ki{a2rip2 + 01(02 — d2)) > gb{Ki) = aiKi{a2 — d2) > 0. 
Therefore, we have gb{xi) > 0 for all xi G [ 0 , Ki] if the following conditions hold 

4AriOi(o2 - oi)(d2 - 02) 

d2 < -;- 

ri (Ari02 — ATiOi -I- 02) 

The discussion so far also indicates that we have both fb{x2) > 0 for X2 G [0,^12] and gb{xi) > 0 for 

xi G [ 0 , ATi] if the following inequalities hold for z = 1 , j = 2 or z = 2,j = 1 


r j 1 ^ 4^1^02(0^ 

Oi > maxjoj, dij, pj < ' - 

n {K^ttj 


ai){dj - Qj) 
KiU^ -I- Oj)^ 
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Now assume that these conditions hold, then we have F{x2) and G(xi) are positive on their restricted 
domain. By algebraic calculations, if Ui > max{(ii,d2} for both i = 1 , 2 , then both F(x2) and G'(xi) have 
its unique critical points xf,i = 1,2 in their restricted domain where 


(nP2 + ai -d2 - 1/(01 - d2)(rip2 + ai -^ 2 )) 

e(o,ifi) 


G ( 0 ,K 2 ). 

If F{x 2) < Ki and G'(xi) < K2, then we can conclude that both maps xi = F{x2) and X2 = G{xi) are 
unimode and the skew product of F x G maps [ 0 , K2] x [ 0 , Ki] to its compact subset. Since both F and G are 
continuous and differentiable, therefore, Xi = F(x2) and X2 = G(xi) has at least one positive intersection for 
X2 G [0,1^2]) G Now we focus on sufficient condition that leads to < Ki and G{x\) < K2 

when oi >02. Since 


riP 2 


and 


K2 {r2Pi + 02 - di - 1/(02 - di){r2Pi + 02 - 


max 

Q<X 2 <K 2 


{F(x 2 )} = F(x^) < 


maxo<a;^<K 2 {Mx 2 )} 

xa.mo<x 2 <K 2 {fb{x 2 )} 


fb{K 2 / 2 ) 

fb{K 2 ) 


K 2 r 2 Pi + 4 di 
4 (oi - di) 


and 


^ xaaxo<a:i<Ki{gt{xi)} 

max {G(xi)} = G(xi) < — - ^ 

0 <Xi<ivi 


5 t(Fi/ 2 ) 


aiFi(riiGip 2/4 + ^2) 


riP 2 (Kiai—JCia 2 —a 2 )^ ’ 


mino<^i<K,{< 7 o(xi)} K,a,{a 2 - d^) - 

therefore, we have F{x2) < Ki and G(xi) < K2 when oi > 02 if the following inequalities hold 

K 2 r 2 Pi + 4 di 4 (F:iOi - Kidi - di) 

< Ki Pi < - 


and 


4 (oi - di) 
K2 > 
P 2 < 


K 2 r 2 


aiKi{riKip 2 / 4 :+d 2 ) 


4 :Kiai{K 2 a '2 — K 2 d 2 —d 2 ) 


airiK^-\-r2K2{Kiai — Kia2—a2)^ ' 

Therefore, we can conclude that Model ( 1 ) has at least one interior equilibrium (x*, yl, x/ yj) if the following 
inequalities hold 

4 {Kiai — Kidi — di) 


li > max{oj, di, ^2}, Oj > max{di,d2}, pi < 


^j '^3 


and 


. r AKiai{aj — ai)(dj — aj) AKiai{Kjaj — Kjdj — dj) 

Pj < mm I ^ ^ ’ aiuKf + r^FT/FTiOi - - 0/2 




In addition, since both F{x2) and G(xi) are unimode maps in their domain with unique local maximum, 
thus, we have 


a:i = F(x 2 ) > F( 0 ) = 


02(^1 


Oir2Pl + O1O2 — 02^1 


and X2 = G(xi) > G( 0 ) = 


aid 2 


O2X1P2 + 0102 — 01^2 


Therefore, we have 


aidj 


jFipj -\-aiaj —aidj 


< X* < Kj for both i = I,j = 2 and i = 2 ,j = 1 . 
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Now assume that oi = 02 = a, then both /& and gb are reduced to linear decreasing functions, i.e., 

= a [r2Pi{K2 - X2) + K2{a - di)] and = a [riP2{Ki - xi) + Ki{a - dz)] 

which indicates that fb{x2) < 0 for X2 G [ 0 ,^ 412 ] if di > a + r2pi and gb{xi) < 0 for xi S [ 0 , ifi] if 
d2 > a + rip2. Therefore, if ai = 02 = a and either di > a + r2Pi or d2 > a + rip2 holds, then Model 
( 1 ) has no interior equilibrium. On the other hand, both fb{x2) > 0 for X2 £ [ 0 ,df 2 ] and gb(xi) > 0 for 
Xi G [ 0 , Ki] if a > max{di, d2}. Then apply the discussions for the case ai ^ 02, then we can conclude that 
Model ( 1 ) has at least one interior equilibrium if 

a > max{di,d2}, Fi^x^) < Ki, and G{xl) < K2. 


Applying the similar arguments for the case Oi > aj, we have 


max {J^(a;2)} = F{x2) < 

0 <X 2 <K 2 


uiaxo<^ 2 <K 2 {ft{x 2 )} 

mino<x 2 <if 2 {/fe( 2 ; 2 )} 


MK2/2) 

fb{K 2 ) 


K 2 '<’ 2 Pi + 4 di 
4 (ai - di) 


and 

max {G(xO} = Gixl) < ^ ^ K^np^ + 4 ^ 

0<xi<Ki mmo<j:^<Ki{9b(xi)} 9b (Fi) 4(a2 - d2) 

Therefore, we can conclude that Model ( 1 ) has at least one interior equilibrium (x*,y^, a;2i 2/2) if fii® following 
inequalities hold 


ai = 02 = a > maxjdi, d2}, Pi < 


4 (A:ia - ATjdi - di) 




for both i = l,j = 2 and i = 2 ,j = 1 . 
^ = 2 J = 1 . 


In addition, 


< X* < Kj hold for both i = l,j = 2 and 


□ 


Proof of Theorem 3.6 

Proof Suppose that min{ai,a2} > d, then we have pi = is the 

unique interior equilibrium of the single patch model ( 2 ) in the absence of the dispersal in predator for both 
1 = 1 , 2 . Now recall from the null dines ( 5 ), we have 


X,{Xj) 


_ PiljGpPjGP+dj _ _ Piqjixj)pj(xj)+di 

ai{l+Piqj{xj))-(piqj(xj)pj{xj)+di) piqj(xj)[ai-pj{xj)]+ai-di 


which indicates that 

^ Piqj{pj)[ai - pj{pj)] + Oi - d Oi-d 

This implies that Xi = pi for i = 1,2 is the positive solutions of the null dines ( 5 ). Therefore, we can solve 
that (pi, ^1,/i2) ^^2) is an interior solution of Model ( 1 ) if min{ai,a2} > d = di = d2. By substituting the 
equilibrium (pi, p2) *^2) into the Jacobian matrix ( 10 ), we obtain its characteristic equation as follows: 

Jd(A) = A"^ + (oi + q; 2 )A^ + [q!iq ;2 + d{f 3 \ + / 32 )]A^ + d{a\j 32 + q; 2 / 3 i)A + d^(/ 3 i /32 ~ 7172 = 0 


where 


(Xi 

r,pPK,a,-K^d-a^-d) ^ > Q K.Oi K,d tt, d < 0 < p, < K, 

li 

_ pii 2 iVj{aj-df ^ Q 

P1P2 - 7172 

1^11^2(0.1—d)^a2 — d)^ {i^iP2-\-i^2pi-\-^) ^ n 
aia 2 
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Then the real parts of the solutions of H(X) are all negative if oi + 02 > 0 while the solutions of H{X) has 
positive if ai + 02 < 0 . Notice that the single patch i ( 2 ) has global stability at (/ii, t'i) if Oi > 0 < 

fj,i < Ki. Therefore, the interior equilibrium is locally asymptotically stable of < Ki for both 

1 = 1,2 while it is unstable if 


— Kid — Old) 
Kiai 


rn 2 {K 2 a 2 - K 2 d - 02 - d) ^ ^ 

1^202 


Assume that 0102 < 0 and Oi + 02 > 0 . Then the real parts of all solutions of id(A) can be still negative 
if 0102 + d(/ 3 i + (32) > 0 and d(ai /32 + 02 / 3 i). By algebraic calculations, we can conclude that if oi < 0 and 
the dispersal of predator yi is large enough, then the interior equilibrium (^1, i^i,/i2,2^2) can still be locally 
asymptotically stable, where pi should satisfy the following condition: 


Pi > max ■ 


—Vj — rj pipj{Kiai — Kid — ai — d)(Kjaj — Kjd — aj — d) 
{KiKjajVjdviiai — d)'^) 


liiUjKj{uiPj-\-\')(^aj—d) [Kiai—Kid — ai—d) ^ 


}■ 


Now if oi = 02 = o, ri = r2 = 1 , ATi = K2 = AT, di = d2 = d. The discussions above implies that Model 
( 1 ) has the same stability at A® = {pL^v,p,v) as the stability of the single patch model ( 2 ) at (/i, i^) where 
p = and V = ■ Therefore, A® is locally asymptotically stable if < p < K while it is 

unstable \i p> 

Now we should show that Model (1) has the unique A® = [p, v, p, v) whenever a > d. Notice that F{x2) 
and G{xi) have the following properties in the symmetric case (i.e., Oi = 02 = a,ri = r2 = i,Ki = K2 = 
K, di = d 2 = d): 

1 . F{x 2) = ppK-x2)+aK(i-d) with F{ 0 ) = and F{p) = F{K) = p. In addition, ^(0:2) > 0 for all 

X2 € [ 0 , AT] and F(x2) has a unique critical point X2 for X2 € [ 0 , AT] where 


P<X 2 


K {^pi + a — d— 1/(0 — d)(pi + o — d)^ 


Pi 


e (o,if). 


2. G{xi) = with G( 0 ) = and G{p) = G(K) = p. In addition, G{xi) > 0 for all 

Xi S [ 0 , AT] and G{xi) has a unique critical point xl for xi S [ 0 , AT] where 

AT (p2 + a — d — {a — d){p2 + a — d )) 

p < xl = - — e ( 0 , AT). 

P2 


The discussions above indicate that both F{x2) and G{xi) are unimode maps with a unique interception at 
Xi = X 2 = p- 

□ 


Proof of Theorem 4-1 

Proof Proof of Item 1 can be obtained by adopting the proof provided in Theorem 2 . 1 . We omit details. 

The stability of Agoooi ^XiOOOj ^00^201 Fk^qk^o can be obtained from eigenvalues of the Jacobian matrix 
of Model ( 7 ) evaluated at these equilibria through simple algebraic calculations. We omit details. But we 
will return to the local stability of when we prove Item 4 . 

Item 3 (a): If Xi = 0 , then we have limt_>oo di(i) = limt_,.oo y2{t) = 0 . This implies that the omega limit 
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set of Model (8) is yi = y 2 = 0. Therefore, prey Xi persists by applying Theorem 2.5 of [32] since 


dxi 


Xidt 


Xi—0 


= ri> 0. 


Item 3(b): Recall Pi{x) = and qi{x) = fQj. j — 1^2. Then we construct the following 

Lyapunov function 


„.(C) — n (K ) 

Viji^i-jUi^yj) = {pj + dj) I T-T d^ + {pj + dj)yi + Piyj- 

JKi Pi\^) 

If Pi > Ki, then we have 

— - = {Pj + dj) [ipiix,) - PiiK,)] qi{xi) + yi{pi{Ki) - d,) < 0 

since Pi{Ki) — di < 0 ^ pi > Ki. Therefore, Model (8) has global stability at {Ki, 0, 0) if pi > Ki. 


Item 3(c): We construct the following Lyapunov function 

P^{0 -P^{P^) 


Vij {xi ,yi,yj) — {pj + dj ) 


PiiO 


+ iPj + dj) 


Vi - Vj 
Vz 


dvi + Pi 


rvj ijj - 

£>• dj 


dvj. 


If 


Ki-l 


< Pi < Ki, then we have 


dt 


{pj + dj)yiyj 


Therefore, Model (8) has global stability at {pi,Vi,v'j) if ^ < pi < Ki. 


Item 4: If Xj = 0, then Model (7) reduces to Model (8) who has global stability at {Ki, 0, 0) if pi > Ki 
while has {pi, Pi, Pj) if < pi < Ki. Therefore, by applying Theorem 2.5 of [32], prey Xj persists if 


dxj 

Xjdt 


Xj—0,yj—0 


= rj > 0 when pi > Ki 


and 


dxi 


Xjdt 


= Tj — QjVj > 0 when 


1 


< Pi < Ki. 


^j=o,yj=i>'. - - 2 

The persistence of both prey can be easily obtained from the persistence of one prey. 


If 


Ki-l 


< Pi < Ki and rj — ajv) < 0, then 


3^j 


Xj=0,Vj=i>' 


= rj — ajpj < 0. According to Theorem 2.18 


by [33], we can conclude that the boundary equilibrium {pi,Vi,0,Vj) is locally asymptotically stable. This 
proves the stability condition of A]' of Item 2. 


Item 5: We construct the following Lyapunov function 

PjiQ-PjiN) 


fVi 


V{xi,yi,Xj,yj) = {pi + di) 


III j 


Pjiij) 


d^j + {pi + di) 


dj 


dj 


-dvj + p^Xi + 


rvz _ pJ 


dVi- 
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Then we have 

^ [{Pi{xi) - Pi(AO] [<li{xi) - H - + PjPiixi) [qi{xi) - vf] 


(Pj+dj)yiyj 

< {pj + dj) [(Piixi) - Piifii)] Mxi) - Vi] - + PjPiixi) - vl] . 

PiPi((pj+dj)yj-Pjyif I ^ „ 

{Pj+dj)yiyj + PjP 


= {pj + dj) [{piixi) - Piifii)] [qiixi) - Vi] - 


- < 0. This implies that 


Therefore, if ^— < flj < Kj and ’’* 4 ^*/^. < hold, then we have - 

Model ( 8 ) has global stability at {0,vil, f2j,Vj). 

Item 6 : Define Viyi,y 2 ) = P 2 yi + Piy 2 , then we have 

^ = P2{piixi) - di)yi + pi{p2{x2) - d2)y2- 

Notice that limsupt_^go a:i(t) < Ki for both i = 1,2. Then if pi > Ki for both i = 1,2, then we have 
max{pi(i4ri) — di,p 2 {K 2 ) — ^ 2 } < —<5 < 0. This implies that 

^ = P 2 iPiixi) - di)yi + pi(p 2 (a: 2 ) - ^ 2)^2 < -Sip 2 yi + Piy 2 )- 

Therefore, both predators go extinct if pi > Ki for both i = 1 and i = 2. Since both limsup(_,,(^ Viit) = 0 
for both i = 1, 2. Then we have Model (7) reduced to the following uncoupled prey model 


x'i = TiXi ( 1 - 


K, 


which converges to Xi = Ki. Thus, Model (8) has global stability at (1411,0,1^2,0) when pi > Ki for both 

i = 1,2. 

On the other hand, ii Q < pi < Kj for both * = 1,2, then we have max{pi(l4ri) — di,p 2 {K 2 ) — ^ 2 } > <5 > 0. 
This implies that 

-jr = P2{piixi) - di)yi + piip2ix2) - ^ 2 ) 2/2 > Sip2yi + Piy2)- 
dt 

Therefore, both predators persist if 0 < < IsTi for both i = 1 and i = 2. 


Item 7 can be obtained from Item 4 and Item 6. 

Item 8 can be obtained from eigenvalues of the Jacobian matrix of Model (7) evaluated at the symmetric 
interior equilibrium {p, v, p, v) through simple algebraic calculations. We omit details. The global stability 
of ip, V, p, v) when < p < K can be obtained by constructing the following Lyapunov function 


Vixi,yi,X2,y2) =P2 

U 

which gives 


Pi(6) -Pi(m) 

p Pi (6) 


P2ii2)-P2 ip) , 

P2(6) 


PIC I I P 2 i^ 2 ) - P 2 ip) , r^V 2 -l 2 

d^i + P2 J — -d?7i + Pi J - s, - d^2 + pi I — ;;- dri2 


ry 
J V 




dVixi Xo 'Uq) PlP 2 ^{.yi ^2) ( 112/ 

’ - = P2iPli^l)-plip))iqiixi)-K)+piip2i^2)-P2ip))iq2ix2)-K)-\ - 


dt 


ym 


□ 
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